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Abstract. We examine the dynamics of a bright solitary wave in the presence of a repulsive or attractive 
localized "impurity" in Bose- Einstein condensates (BECs). We study the generation and stability of a 
pair of steady states in the vicinity of the impurity as the impurity strength is varied. These two new 
steady states, one stable and one unstable, disappear through a saddle-node bifurcation as the strength of 
the impurity is decreased. The dynamics of the soliton is also examined in all the cases (including cases 
where the soliton is offset from one of the relevant fixed points). The numerical results are corroborated 
by theoretical calculations which are in very good agreement with the numerical findings. 

PACS. 03.75,-b Matter waves - 52.35.Mw Nonlinear phenomena: waves, wave propagation, and other 
interactions 



1 Introduction 

In the past few years, the rapid experimental and theoreti- 
cal developments in the field of Bose-Einstein condensates 
(BECs) pQ have led to a surge of interest in the study 
of the nonlinear matter waves that appear in this con- 
text. More specifically, experiments have yielded bright 
solitons in self-attractive condensates ( 7 Li) in a nearly 
one-dimensional setting [2], as well as their dark |3] and, 
more recently, gap 0] counterparts in repulsive conden- 
sates, such as 87 Rb. The study of these matter-wave soli- 
tons, apart from being a topic of interest in its own right, 
may also have important applications. For instance, a soli- 
ton may be transferred and manipulated similarly to what 
has been recently shown, experimentally and theoretically, 
for BECs in magnetic waveguides [Hj and atom chips pj. 
Furthermore, similarities between matter and light waves 
suggest that some of the technology developed for optical 
solitons may be adjusted for manipulations with MWs, 
and thus applied to the rapidly evolving field of quantum 
atom optics (see, e.g., [H]). 

One of the topics of interest in this context is how 
matter- waves can be steered/manipulated by means of 
external potentials, currently available experimentally. In 

a http://nlds.sdsu.edu/ 



addition to the commonly known magnetic trapping of 
the atoms in a parabolic potential, it is also experimen- 
tally feasible to have a sharply focused laser beam, such as 
ones already used to engineer desired density distributions 
of BECs in experiments Depending on whether it is 
blue-detuned or red-detuned, this beam repels or attracts 
atoms, thus generating a localized "defect" that can in- 
duce various types of the interaction with solitary matter 
waves. This possibility was developed to some extent in 
theoretical [H] and experimental studies of dynamical 
effects produced by moving defects, such as the generation 
of gray solitons and sound waves in one dimension 
and formation of vortices in two dimensions (see, e.g., 
and references therein). The interaction of dark solitons 
with a localized impurity was also studied l.'i . 

Our aim here is to examine the interaction of a bright 
matter- wave soliton with a strongly localized (in fact, 5- 
like) defect, in the presence of the magnetic trap. Our 
approach is different from that of Ref. ^Sj , m that we will 
view the presence of the defect as a bifurcation problem. 
We demonstrate that the localized perturbation (indepen- 
dently of whether it is attractive or repulsive) creates an 
effective potential that results in two additional localized 
states (one of which is naturally stable, while the other is 
always unstable) for sufficiently large impurity strength. 
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As one may expect on grounds of the general bifurcation 
theory, these states will disappear, "annihilating" with 
each other, as the strength of the impurity is decreased 
below a threshold value. We will describe this saddle-node 
bifurcation in the present context. We will also compare 
our numerical predictions for its occurrence with analyti- 
cal results following from an approximation that treats the 
soliton as a quasi-particle moving an effective potential. 
Very good agreement between the analytical and numeri- 
cal results will be demonstrated. Finally, we will examine 
the dynamics of solitons inside the combined potential, 
jointly created by the magnetic trap and the localized de- 
fect. Both equilibrium positions and motion of the free 
soliton will be considered in the latter case. 

The paper is structured as follows: in Section we 
present our effective potential theory. In Section we 
discuss numerical methods and results, and provide their 
comparison with the analytical predictions. Finally, in Sec- 
tion m we summarize our findings and present our conclu- 
sions. 



2 Setup and Theoretical Results 

In the mean-field approximation, the single-atom wave- 
function for a dilute gas of ultra-cold atoms very accu- 
rately obeys the Gross-Pitaevskii equation (GPE). Al- 
though the GPE naturally arises in the three-dimensional 
(3D) settings, it has been shown |14U15lllfi| that it can be 
reduced to its one-dimensional (ID) counterpart for the 
so-called cigar-shaped condensates. Cigar-shaped conden- 
sates are created when two transverse directions of the 
atomic cloud arc tightly confined, and the condensate is 
effectively rendered one-dimensional, by suppressing dy- 
namics in the transverse directions. The effective equa- 
tion describing this quasi-lD is simply tantamount to a 
directly written ID GPE. In normalized units, it takes 
the well-known from, 

iiM = -^u xx + g\u\ 2 u + V(x)u, (1) 

where subscripts denote partial derivatives and u(x,t) is 
the one-dimensional mean-field wave function. The nor- 
malized ID atomic density is given by n = \u(x, t)\ 2 , while 
the total number of atoms is proportional to the norm of 
the normalized wavefunction u(x,t), which is an integral 
of motion of Eq. JU : 

/+oo 
\u(x,t)\ 2 dx. (2) 
-oo 

The nonlinear coefficient in Eq. l(T)l is a = zbl. for repulsive 
or attractive interatomic interactions respectively. Finally, 
the magnetic trap, together with the localized defect, are 
described by a combined potential V(x) of the form 

v(x) = ~n 2 x 2 - V S(x - o. (3) 

In Eqs. Q and Q, the space variable x is given in 
units of the healing length £ = h/ ^/noguyfn, where no 



is the peak density, and the normalized atomic density 
is measured in units of no. Here, the nonlinear coeffi- 
cient is considered to have an effectively ID form, namely 
gio = 53£)/(27rZj_), where g^D — 4:7rh 2 a/m is the original 
3D interaction strength (a is the scattering length, m is 
the atomic mass, and l±_ = \Jn~Jmjjj\_ is the transverse 
harmonic-oscillator length, with to± being the transverse- 
confinement frequency). Further, time t is given in units 
of £/c (where c = \J n$gir> Jm is the Bogoliubov speed of 
sound), and the energy is measured in units of the chem- 
ical potential, /i = <?id^o- Accordingly, the dimensionless 
parameter J? = fu^x/giDfio (where oj x is the confining 
frequency in the axial direction) determines the effective 
strength of the magnetic trap in the ID rescaled equations. 
Positive and negative values of Vq corresponds, respec- 
tively, to the attractive and repulsive defect. Finally, since 
we are interested in bright matter- wave solitons, which 
exist in the case of attraction, we hereafter set the nor- 
malized nonlinear coefficient g = — 1. 

It is worth mentioning that modified versions of the 
ID GPE are known too. One of them features a non- 
polynomial nonlinearity, instead of the cubic one in Eq. 
||TJ 16 . A different equation was derived for a case of 
a very strong nonlinearity, so that the local value of the 
potential energy exceeds the transverse kinetic energy. It 
amounts to the same cubic equation ljT)l. but with a non- 
canonical normalization condition, with the integral in Eq. 
© replaced by f*™ \u(x, i)| 4 dx. 

In the absence of potential, Eq. supports stationary 
soliton solutions of the form 

u s (x) = n sech [t](x — £)] exp (in 2 i/2) (4) 

where n is an arbitrary amplitude and £ is the position 
of the soliton's center. It is possible to generate mov- 
ing solitons (with constant velocity) by application of the 
Galilean boost to the stationary soliton in Eq. J2J. 

One can examine the persistence and dynamics of the 
bright solitary waves in the presence of the potential V(x) 
by means of the standard perturbation theory (see e.g., 
Ref. |T71T8| and a more rigorous approach, based on the 
Lyapunov-Schmidt reduction, that was developed in Ref. 
|19p. This method, which treats the soliton as a particle, 
yields effective potential forces acting on the particle from 
the defect and the magnetic trap, 

F dcf = 2n 3 tanh(n(e - C))scch 2 (n(£ - Q)V , (5) 

F trap = -2r? 2 Cn, (6) 
which enter the equation of motion for £(i): 

C -^imp ~t~ -ftmp. (7) 

Below, results following from this equation will be com- 
pared to direct simulations of Eq. . 

The stationary version of Eq. J7J) (£ = 0) , 

(n 3 Vb/fi 2 ) (tanh 9) scch 2 6» = ^ - 6, (8) 
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with 9 = r](£ — £), determines equilibrium positions (£) of 
the soliton's center. Depending on parameters, this equa- 
tion may have one or three physical solutions, see below. 
In what follows, we will examine the solutions in detail 
and compare them to numerical results stemming from 
direct simulations of Eq. . 



3 Numerical Methods and Results 

In order to numerically identify standing wave solutions of 
Eq. QJ, we substitute, as usual, u(x,t) — exp(iAt)w(x), 
which results in the steady-state problem: 



Aw = ^w xx 



w 



V(x)w. 



(9) 



This equation is solved by a fixed-point iterative scheme 
on a fine finite-difference grid. Then, we analyze the stabil- 
ity of the obtained solutions by using the following ansatz 
for the perturbation, 



u(x) 



JAt 



w(x) + a(x)e 



-At 



b*{x)e 



-\*t 



(10) 



(the asterisk stands for the complex conjugation), and 
solving the resulting linearized equations for the pertur- 
bation eigenmodes {a(x),b(x)} and eigenvalues A associ- 
ated with them. The resulting solutions are also used to 
construct initial conditions for direct numerical simula- 
tions of Eq. (JIJ, to examine typical scenarios of the full 
dynamical evolution. To eliminate effects of the radiation 
backscattering in these simulations, we have used absorb- 
ing boundary layers, by adding a term in Eq. (Ip. of the 
form: 

r u = - [(1 + tanh(1000(a; - R)) 

(11) 

+ (1 - tanh(1000(x - L))} u, 

which is defined on the domain L < x < R. In the nu- 
merical simulations presented herein, the (^-function of the 
potential was approximated by a Gaussian waveform, ac- 
cording to the well-known formula, 



5(x) 



lim 

T^0 + 



1 



2ira 



ex P -T- 



4a 



(12) 



Lorentzian and hyperbolic-function approximations to the 
(5-function were also used, without producing any conspic- 
uous difference in the results. 

As mentioned above, depending on the value of the 
defect's strength Vq, Eq. JHJ may have either one or three 
physical roots for £ (the equilibrium position of the soli- 
ton's center) . The border between these two generic cases 
is a separatrix where two of the roots merge in one be- 
fore they disappear. All the qualitatively different cases 
are depicted in the top-left panel of Fig. ^ The physical 
interpretation of this result can be given as follows. Ob- 
viously, in the absence of the defect there exists a stable 
solitary- wave configuration centered at C = (hence there 
is a single steady state in the problem). On the other hand, 
it is easy to see that Eq. (JSJ) generates three solutions for 
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Fig. 1. Saddle-node bifurcation of stationary states for the lo- 
cation £ of the bright soliton inside the magnetic trap (S7 = 
0.1), with the localized defect of strength Vo located at £ = 6. 
The top-left panel displays the corresponding solutions of the 
stationary equation JHJ. For the weak defect, dash-dotted line, 
only one steady state exists very close to the origin. As the 
defect's strength increases, two additional fixed points (both 
located on the same side of the impurity) are created in a 
saddle-node bifurcation. The top-right panel depicts the posi- 
tion and stability (solid for stable and dashed for unstable) for 
the steady states as a function of the defect's strength Vo- The 
thin horizontal line for £ = 6 shows the location of the defect. 
The two bottom plots depict another version of the stability 
diagram, in terms of the soliton's norm P [see Eq. ©], as Vo 
is varied. 



large Vo- Hence, there should be a bifurcation point, of the 
saddle-node type, that leads to the disappearance of two 
branches of the solutions as Vo decreases. Furthermore, 
based on general bifurcation theory principles, one of the 
steady states disappearing as a result of the bifurcation 
may correspond to a stable soliton, whereas its compan- 
ion branch definitely represents an unstable solitary wave. 
The full bifurcation-diagram scenario, for the position of 
the soliton's center and its norm, is depicted in Fig. ^ It 
is interesting to note that the norm of the solitons corre- 
sponding to the unstable branches varies almost linearly 
with the defect's strength, while the stable branches corre- 
spond to solitons whose norm is approximately constant. 

The qualitative predictions about the nature of the 
steady states and their stability have been tested for re- 
pulsive and attractive defects, as shown in Figs. [5] and [2 
respectively. In these figures, three left panels show the 
spatial profiles of the stable branch at ( = 0, and the 
unstable and stable branches in the neighborhood of the 
defect. The middle panels show the temporal evolution of 
each one of these solutions, while the right panels show 
the results of the linear stability analysis. The latter set 
clearly illustrates the instability of the middle branch due 
to the presence of a real eigenvalue pair. It is also note- 
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worthy that, in the case of the repulsive defect (in which 
case an unstable solution is centered at the defect) the 
soliton oscillates around the nearby stable steady state, 
shedding radiation waves, cf. Fig. [5] On the other hand, 
in the attractive case the unstable solution centered be- 
side the defect is "captured" by the defect, resulting in 
its trapping at the defect's center, cf. Fig. [3] However, a 
fraction of the condensate is also emitted from the defect 
in the process, leading to oscillations that can be observed 
in the respective space-time-evolution panel. 

To verify the analytical results following from Eqs. JJJ) 
and JHJ), we have compared the analytically predicted crit- 
ical value of Vo (for which a double root appears) with the 
numerically obtained turning point for the saddle-node bi- 
furcation. This comparison was performed for many val- 
ues of the impurity center £. In fact, the critical value was 
predicted using two different forms of the analytical pre- 
diction: one with the Dirac <5-function proper, and another 
one with the Gaussian approximation for the (^-function 
and an accordingly modified version of Eq. JSJl, namely 

oo 

F ilap +F tTap = J V(x)-^(\u(x)f}dx = 0, (13) 

— oo 

with the function V{x) incorporating the parabolic mag- 
netic trap and the Gaussian impurity terms. Here the inte- 
gration was performed with the numerically implemented 
V(x), and the best fit of u(x) to a hyperbolic secant wave- 
form has been used in Eq. i|13|) . The parameter values 
along with the resulting critical values of Vq are given in 
Fig. 0] In all cases, the numerical results for the bifurca- 
tion point closely match the theoretical predictions. 

Having examined statics and dynamics in the vicinity 
of the stable and unstable fixed points of the system, we 
now turn to an investigation of the dynamics, setting the 
initial soliton farther away from the equilibrium positions. 
Figure displays three typical examples, with the soliton 
set to the left and to the right of the repulsive (and of the 
attractive) defect. In the repulsive case, we observe that 
the soliton is primarily reflected from the defect; however, 
when it has large kinetic energy at impact (which takes 
place if it was initially put at a position with large po- 
tential energy), a fraction of the matter is transmitted 
through the defect. On the other hand, in the attractive 
case, a fraction of the matter is always trapped by the 
defect. However, this fraction is smaller when the kinetic 
energy at impact is larger. 

We note in passing that, while Eq. Q) can predict not 
only the equilibrium positions of the soliton but also the 
dynamical behavior of £(i), we have opted not to use it 
for the numerical experiments. The main reason is that, 
as can be clearly inferred from Figs. 12131 and 151 the inter- 
action of the solitary wave with the defect entails emis- 
sion of a sizable fraction of matter in the form of small- 
amplitude waves, which, in turn, may interfere with the 
solitary wave and significantly alter his motion (see e.g., 
Fig. |5J). Hence, the prediction of the dynamics based on 
the adiabatic approximation, which is implied in Eq. (|7|). 
would be inadequate in the presence of these phenomena. 



4 Conclusions 

In this work, we have examined the interaction of bright 
solitary waves with localized defects in the presence of 
magnetic trapping, which is relevant to Bose-Einstein con- 
densates with negative scattering lengths. We have found 
that the defect induces, if its strength is sufficiently large, 
the existence of two additional steady states (bifurcating 
into existence through a saddle-node bifurcation), one of 
which is stable and one unstable. We have constructed 
the relevant bifurcation diagram and explicitly found both 
the stable and the unstable solutions, and quantified the 
instability of the latter via the presence of a real eigen- 
value pair. The dynamical instability of these unstable 
states leads to oscillations around (for repulsive defects) 
and/or trapping at (for attractive defects) the nearby sta- 
ble steady state. Additionally, we have developed a collective- 
coordinate approximation to explain the steady soliton 
solutions and the corresponding bifurcation. We have il- 
lustrated the numerical accuracy of the analytical approx- 
imation by comparison with direct numerical results. We 
have also displayed, through direct numerical simulations, 
the dynamics which follows setting the initial soliton off a 
steady-state location. Noteworthy phenomena that occur 
in this case are the emission of radiation by the soliton 
colliding with the repulsive defect, and capture of a part 
of the matter by the attractive one. 

These results may be relevant to the trapping, manip- 
ulation and guidance of solitary waves in the context of 
BEC. They illustrate the potential of the combined effect 
of magnetic and optical (provided by a focused laser beam) 
trapping to capture (either at or near the laser-beam- 
induced local defect) a solitary wave which can be subse- 
quently guided, essentially at will. Naturally, the beam's 
intensity must exceed a critical value, which can be explic- 
itly calculated in the framework of the developed theory. It 
would be particularly interesting to examine the predicted 
soliton dynamics in BEC experiments. 
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Fig. 2. Steady states of the bright soliton for the repulsive defect: Vb = — 1, a — 0.045, r\ = \[2, fi = 0.1 and £ = 6. The first 
row corresponds to the steady state at £ = 0, the second row to the steady state centered at the defect, and the third row to 
the steady state trapped to the right of the impurity. For each row, the left panel displays the numerically exact steady-state 
soliton profile, the middle panel is the space-time evolution shown by means of contour plots, and the right graph shows the 
spectral plane (A r , Ai) for the stability eigenvalues A = A r + iA, corresponding to this solution. For the two stable steady states 
(trapped at the defect and to the right of it), the solution remains stationary as expected. On the other hand, for the unstable 
steady state, after approximately 20 time units, the instability fragments the soliton into a more localized part oscillating to its 
right and a more extended part oscillating to its left. 
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Fig. 4. Critical values of Vb, corresponding to the disappearance of two steady states, for a = 0.045. The left and right graphs 
pertain to the attractive and repulsive defect, respectively. 
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Fig. 5. Examples of the soliton interaction with the repulsive defect (top panels) and the attractive one (bottom panels) located 
at £ = 6. The soliton is initially offset with respect to the steady states of the model. For all cases, the parameters are Vo = ±1, 
a = 0.045 and r\ = v2. The initial position of the soliton is £0 = 12, —6, —12, for the plots from left to right. In the repulsive case 
(top panels), the soliton is primarily reflected by the defect (with a small transmitted fraction of the norm). Similar behavior 
was observed for other values of f, with the amount of material passing through the defect increasing with £o- On the other 
hand, in the attractive case (bottom panels), the soliton gets fragmented into reflected, trapped and transmitted parts. For 
larger initial values of f , the trapped fraction is smaller. 



